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FOREWORD 


This research report is concerned with the dynamic response of unidirectional 
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C. Chamis who has carefully reviewed this report and provided a number of concrete 
suggestions. 
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OFF-AXIS IMPACT OF UNIDIRECTIONAL COMPOSITES WITH CRACKS: 
DYNAMIC STRESS INTENSIFICATION 


by 

G. C. Sih 

Institute of Fracture and Solid Mechanics 
Lehigh University 
Bethlehem, Pennsylvania 18015 

and 

•k 

E. P. Chen 
Sandia Laboratories 
Albuquerque, NewMexico 87115 

ABSTRACT 

The dynamic response of unidirectional composites under off-axis (angle 
loading) impact is analyzed by assuming that the composite contains an initial 
flaw in the matrix material. Because of the complexities that arise from the 
interaction of waves scattered by the crack with those reflected by the inter- 
faces within the composite, dynamic analyses of composites with cracks have been 
treated only for a few simple cases. One of the objectives of the present work 
is to develop an effective analytical method for determining dynamic stress so- 
lutions. This will not only lead to an in-depth understanding of the failure 
of composites due to impact but also provide reliable solutions that can guide 
the development of numerical methods. 

The analysis method utilizes Fourier transform for the space variable and 
Laplace transform for the time variable. The time-dependent angle loading is 


*This work was completed during Dr. Chen's tenure at Lehigh University. 
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separated into two parts: one being symmetric and the other skew- symmetric with 

reference to the crack plane. By means of superposition, the transient boundary 
conditions consist of applying normal and shear tractions to a crack embedded in 
the matrix of the unidirectional composite. Mathematically, these conditions 
reduce the problem to a system of dual integral equations which are solved in 
the Laplace transform plane for the transform of the dynamic stress intensity 
factor. The time inversion is carried out numerically for various combinations 
of the material properties of the composite and the results are displayed graph- 
ically. 


INTRODUCTION 

Past work on the development of high performance composite materials was 
mainly concerned with achieving high strength and modulus. This requirement 
alone, however, may result in a composite that is excessively brittle and lacks 
the ability to resist impact loading. The energy absorption or toughness of the 
composite is also an important property that must be accounted for in addition 
to strength and stiffness. 

The concept of fracture toughness has mostly been applied to homogeneous 
isotropic materials [1] based on the linear fracture mechanics theories such as 
those advanced by Griffith, Irwin and others. These theories, developed for 
single-phase materials, have had limited success in characterizing the fracture 
behavior of composites which are inherently nonhomogeneous and anisotropic. 

This is mainly because the fracture modes in composites are multi-facet and can 
include interface failure, fiber breaking, matrix fracture, etc. The individ- 
ual contribution of each of these failure modes is not clearly accounted for 
and/or not related to the critical failure load. As a result, large discrepancies 
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between the theory and experiment can result. 


A study on the selection of appropriate mathematical models for different 
unidirectional composite systems was made [2] in the case of static loading. 

Many of the assumptions in [2] will also be used in the dynamic problem treated 
here. One of them is the existence of inherent flaws or cracks which are the 
sites of failure initiation. 

Analytical investigation of the fracture of fibrous composite materials sub- 
jected to impact loading has been meager because the elastodynamic stress analy- 
sis involves numerous parameters and is enormously complex. This is necessitated 
by the complex nature of the dynamic load transfer characteristics in composites 
containing initial imperfections such as flaws or cracks. The stress wave solu- 
tion is not only time-dependent but it interacts with the material properties of 
the constituents of the composite and the various geometric parameters. The in- 
fluence of these parameters will be analyzed in this impact study with particular 
emphasis placed on determining the dynamic stress intensity factors k-j and k2 
arising from normal and shear loading. Their combination (off-axis or angle 
loading) determines the response to loading of a more general nature and reflects 
the energy absorption property of the composite. Several examples of how k-j and 
k 2 can be combined to predict crack behavior in dynamic stress fields are found 
in [3]. The question of whether there is the need of how to define a dynamic 
fracture toughness parameter differing from its corresponding static value has 
been the subject of many past and present debates within the fracture mechanics 
community. Thus far, no general agreement has been achieved. 
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This report is concerned with dynamic fracture analysis and, particularly, 
with the development of an analytical method for obtaining effective dynamic 
stress solutions to unidirectional composites with cracks embedded in the ma- 
trix. Other possible failure modes will be dealt with in future reports. Ef- 
fective stress solutions for and ^2 are essential as they are the prerequi- 
sites for formulating failure criteria and guiding the development of numerical 
procedures. 


ANGLE CRACK UNDER IMPACT 

Figure 1(a) considers a crack in a layer of matrix material of thickness 
2h. The composite is reinforced by unidirectional fibers that are aligned par- 
allel with one another and make an angle with the time-dependent applied stress 
a(t). Without serious loss in generality, the composite is assumed to be modeled 
by a layer of cracked material with elastic properties pp v-j and p-j sandwiched 
in between two dissimilar media with properties ^ 2 * V 2 -and P 2 > Figure 1(b). The 
number of layers surrounding the cracked layer is reasonably large so that the 
average shear modulus pg* Poisson's ratio V 2 and mass density P 2 can be used. 

The basic two-dimensional elastodynamic equations in the theory of elastic- 
ity can be expressed in terms of two scalar potentials (l).(x,y,t) and i|).(x,y,t) 

J J 

where i,j = 1,2 with 1 and 2 referring to the cracked layer and the surrounding 
material, respectively. In terms of the Lamd coefficients X. and p., the dy- 

J J 

namic stress components are 
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"'’ 3x9y^ 

J 


3^<j, d^rli. 

*^0 ^ ^3y^ " 3x3y^ 

J 


X. X.+2y. 

<®z>, = ^ <x^> *j 

J J J 


xy' 


324, 324, 324, 

(t^vx/) ■ (2 0x3y '*’ " 9x^^^ 


(1) 


where = 3^/3x^ + 3^/3y^ and the thickness shear stresses are assumed to van- 
ish. The corresponding in-plane displacements are given by 


(u ) = + !!i 

. 3x 3y 

J 




(2) 


Under plane strain, the material elements are constrained in the z-di recti on. 
The governing differential equations can then be obtained from the equations of 
motion: 



- 1 !!!i 



1 !!!i 


(3) 
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in which c 


Ij 


and 


‘"2j 


are the dilatational and shear wave velocities defined as 


X..+2V. 


1/2 



(4) 


The problem involves the determination of the potentials 4>.(x,y,t) and i|^.(x,y,t) 

J J 

in equations (3) from the transient boundary conditions of the crack problem. 

The analysis may be simplified considerably if the problem is separated in- 
to two parts. The first concerns with normal stresses applied to the crack such, 
that symmetry prevails about the x-axis in Figure 1(b) while the second deals 
with shear surface tractions so that the problem is skew- symmetric with refer- 
ence to the x-axis. Both of these problems will be presented separately. 

NORMAL IMPACT 

Let the composite body be initially at rest such that the stresses are zero 
everywhere. Suddenly, at t=0, a normal stress of magnitude -a^ is applied to 
the top and bottom crack surfaces in Figure 1(b) and kept on the crack of length 
2a thereafter. Referring to the set of axes x and y that are placed parallel 
and normal to the line crack, the following conditions are prescribed: 


(aj (x,o,t) = -o„H(t); (t^,,) (x,o,t) = 0, 0<x<a; t>0 


1 


xy' 


(5) 


1 


where H(t) is the Heaviside unit step function. The symmetry conditions about 
the axis y=0 are enforced by noting 


(u^) (x,o,t).= 0; (t^J (x,o,t) = 0, x^a; t>0 
^.1 1 


( 6 ) 
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Perfect bonding will be assumed along the interfaces between material 1 and ma- 
terial 2. This requires the stresses and displacements to be continuous across 
y = ±h. On account of symmetry, only the upper half plane y^O need to be con- 
sidered, i.e., 

(cTy) (x,h,t) = (a ) (x,h,t) 

y 1 y 2 

(7) 

(t„„) (x,h,t) = Ct^ ) (x,h,t) 
xy 1 xy 2 


and for the stresses and 


(u ) (x,h,t) = (u ) (x,h,t) 

X 1 X 2 

( 8 ) 


(u ) (x,h,t) = (u ) (x,h,t) 

y 1 y 2 

for the displacements. 

Vual ZntzgAxit zqucvUoyu, It is convenient at this point to apply the 
Laplace transform to the time variable t which corresponds to p in the trans- 
formed plane. Consider the standard Laplace transform on f(t): 

f*(p) = / f(t) e'P^dt (9) 


whose inversion is 


= i / f*(p) e'’*dP (10) 
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in which Br stands for the Bromwich path of integration. The application of 
equation (9) to equations (3) yields 



Again, the condition of syiranetry requires only the consideration of x and y in 
the first quadrant. The Fourier cosine and sine transforms defined by 

00 

f^(s) = / f(x) cos(sx)dx 
0 

( 12 ) 

00 

f(x) = f / f^^Cs) cos(sx)ds 
0 

and 


f^(s) = / f(x) sin(sx)dx 
0 

03 ) 

f(x) = |- / f^(s) sin(sx)ds 
0 

are now applied to the space variable x. This simplifies equations (3) to a set 
of ordinary differential equations which can be solved giving 


<|)*(x,y,p) = f / [A^^hs,p)e 

0 


-Yiiy 


.( 2 ) 


Yiiy_ 


+ A' '(s,p)e ] cos(sx)ds 
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(14) 


*f(x,y,p) =1? [B^’Us.P)e + B^^'(s.p)e^^’^3 s1n(sx)ds 


for the cracked matrix and 


4>fCx.y,p) = f / C^^^(s,p)e cos(sx)ds 
0 

^|(s,y,p) = I / C^^^s,p)e sin(sx)ds 


05) 


for the averaged fiber-matrix material. In equations (14) and (15), the quanti- 
ties Y-|j and Y2j are given by 


n2 1/2 n2 1/2 

l^lj ' %> ’ ’^2j = 


(16) 


The functions in equations (14) and (15) are deter- 
mined from the transient boundary conditions. To this end, equations (5) and 
(6) will be written in the Laplace transform plane: 


(cr*) (x,o,p) 
^ 1 


(x*y)^(x,o,p) = 0, 0<x<a 


(17) 


and 


(u*) (x.o,p) = 0; (t* ) (x,o,p) = 0, x>a 


1 


xy' 


( 18 ) 


1 


In the same way, equations (7) become 
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(a*) Cx,h,p) = {a*) tx,h,p) 

y 1 y 2 


(t* ) Cx.h.p) = (t* ) (x,h,p) 
xy ^^2 

and equations (8) take the forms 


09 ) 


(u*) (x,h,p) = (u*) (x,h,p) 

X 1 X 2 

(20) 

(u*) (x,h,p) = (u*) (x,h,p) 

y 1 y 2 

The stresses and displacements in equations (1) and (2) may also be trans- 
formed into the Laplace transform plane. Without going into details, the appro- 
priate Laplace transform of the stress and displacement expressions in equations 
(17) to (20) may be used to satisfy all of the necessary boundary, symmetry and 
continuity conditions. This leads to the following set of dual integral equa- 
tions: 


oo 

/ A(s,p) cos(sx)ds = 0, x^ia 
0 


oo Tta 

/ sFj(s,p) A(s,p) cos(sx)ds = - 


( 21 ) 


in which Fj(s,p) stands for the known function 
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p /c n^ - ^ -rr— r<;^+v2 V ^"'^11^21^ 1 

Fl(s,p) - s(i_K^“ U4 (S +T2V 2 ® ® ** 

^s(s^n|,)[r2,(B'')6(^>-6(^)B<2>)-ri,]e'''^”^2’”’ 

+ tj {s“n|,)' + (22) 

while the quantities k-j and Aj are defined as 

(23) 

„ ^,^,[eC2) , 3 ( 3 ) , -2(^11^21 )\ 3 ( 4 ) / 2 v,,h 

such that 6^'^. — . 8^^ are given by 

gO) = („(3)„(6) . 3(2) . (^(4)„(6) . „(2)„(8))/,^ 

(24) 

3(3) . („(1)„C7) . „(3)„(5))/,^. 3(4) , („(1)„(3) . „(4)„(5))/,^ 
where is 

4„ = u^’)e(«) - e<2)at5> (25) 

The quantities in equations (25) are complicated functions 

of s, p and the material constants. They are given by equations (I.l) in Appendix 


The problem is now reduced to finding the single unknown A(s,p) governed 
by equations (21). Once A(s,p) is known, the functions 
that are required in equations (14) and (15) for the Laplace transform of the 
potentials c|)'^(x,y,p) and 4)^(x,y,p) can be obtained from equations (1.2) outlined 
in Appendix I. What remains is the determination of a solution for the dual in- 
tegral equations (21). This will be accomplished with the help of a method by 
Copson [5] which has been used by Chen and Sih [6] for solving dynamic crack 
problems involving single-phase homogeneous materials. Following the details 
in [5,6], it can be shown that 


A(s,p) = - 


Tra^a 




C26) 


is a solution of equations (21) with being the zero order Bessel function of 
the first kind. The function $|(C,p) is calculated numerically from a Fredholm 
integral equation of the second kind: 

1 

^f(C.p) + / $f(n,p) Kj(S.n,p)dn = (27) 


whose kernel 


Kj(?,n,p) = ^ / s[Fj(f. p) - 1] Jq(sC) Jo(sn)ds (28) 
is symmetric in ? and n* 

Mode. I dynamic. lyitemlly ^actoa. The transmission of the time-de- 
pendent load to the vicinity of the crack tip can be best described by the in- 
tensification of the local stresses. A quantity that has been used widely in 
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the static theory of fracture mechanics is the "stress intensity factor" which 
can be extracted from the asymptotic expansions of the stresses near the crack 
tip. Referring to Figure 2, let r-j and 8^ be a set of local polar coordinates 
measured from the right hand crack tip located at x=a and y=0 in the matrix 
material, the singular character of the dynamic stresses is described only by 
the space variables and hence can be more easily determined in the Laplace 
transform domain. This observation was first made by Sih, Ravera and Embley 
[7]. Following their procedure, the local stresses in terms of r^ and S-] are 
found: 


k?(p) 0, 6-1 30, 

(c^) »P) = '• — cos ^ (1 - sin sin -y) + 0(r| ) 


k*(p) 0, 0, 30, 

(^w) (t'vQvP) = — cos ^ (1 + sin Y- sin + 0(rp 

(29) 

k?Cp) 01 

(o*) (r,,0,,p) = — 2v, cos Y- + O(i^i) 

z 1 I I I 

k*(p) 0, 0, 30, 

CT*y)^(ri,0i,p) = sin ^ cos ^ cos ^ + 0(r|) 

Only the dynamic stress intensity factor, k^(p), in equations (29) need to be 
inverted to real time t: 


k^(p) 


$|(l,p) 

P 



(30) 


where the function $|(l,p) is found from $|(?,p) by letting the nondimensional 
parameter 5=1 representing the crack tip location. The functional dependence 
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of the stresses in r-j and 0^ as shown by equations (29) reveals that the dynamic 
stresses also possess the inverse square root singularity in terms of r-j and 
that the angular distribution in 0, is the same as the case for static loading. 


Applying the Laplace inversion formula in equation (10) to (30) renders 
the factor k^(t) as a function of time, i.e.. 


a Va Of(l,p) t 


(31) 


It is apparent that $|(l,p) must be first known before the integration of equa- 
tion (31) can be performed. Refer to Appendix II for a detailed account of the 
procedure used for evaluating equation (31). Three different sets of $^(l,p) 
values are plotted against the dimensionless Laplace transform wave number C2-|/pa. 
They are given in Figures 3 to 5 for = P2 and v-j = V2 = 0.29 while the ratios 
a/h and p 2 /i^t s'"® varied. In general, all the curves tend to rise quickly and 
then flatten out. It would be more meaningful to discuss the influence of a/h 
and ^2^^! stress intensity factor k-j(t). 

Figures 6 to 8 display the normalized stress intensity factor k^(t)/0Q/a 
as a function of C2-|t/a. In Figure 6, the crack length to layer thickness ratio, 
a/h, is fixed at unity while the shear moduli ratio, y2/y-| ts increased from 0.1 
to 10.0 as indicated. The k^(t) factor is oscillatory in nature reaching a peak 
and then decreases in magnitude as time increases. The oscillation is more pro- 
nounced when the shear modulus of the cracked material is greater than that of 
the surrounding material, i.e., y-i > V2' values of k^(t) decrease below 
those of the corresponding homogeneous case, = P2» solved previously by Chen 
and Sih [6] when < P2. The influence of a/h on k^(t) is exhibited in Figures 

7 and 8 for the two cases of P2/P1 =0.1 and 10.0, respectively. For P2/P1 = 0.1 
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in Figure 7, a decrease in a/h tends to lower the stress intensity factor. Ob- 
served also is a small step in the curve for a/h = 2.0 and small time t. This 
corresponds to the reflection of elastic waves from the material interface. 

The size and time scale are such that this effect showed up quantitatively in 
the graph while the same effect was not noticeable in the other curves. For 
the smaller ratios of a/h such as 0.5 and 1.0, the crack tips are further away 
from the interface and the influence of the reflected waves are not as pro- 
nounced. The opposite trend is observed in Figure 8 for ~ 10.0. When the 

outer material is more rigid than that of the center layer, k-j(t) tends to in- 
crease in magnitude as a/h is decreased. Again, a distinct step in the curve 
for a/h = 2.0 is seen for small time t. As time increases, all of the results 
here reduce to the corresponding static solutions of Hilton and Sih [83. 

GenoAot loading. If the normal stress applied to the crack surface is not 
constant in magnitude but may vary as a function of x, then the dynamic stress 
intensity factor can be obtained by adding a sequence of solutions corresponding 
to step loadings with different stress levels Og, ct-j, etc. In other words, the 
general loading a(t) may be considered as the sum: 

a(t) = cfgHCtg) + cr^H(t^) + a2H(t2) + ... (32) 

This is illustrated graphically in Figure 9. From equations [31) and (32), the 
factor k^(t) that corresponds to a(.t) may be written down immediately as follows: 

, $tn>P) nt 

^ <?-|HCti) + ...] / p e dp (33) 

Equation (.33) may be used to derive k^(t) for any time-dependent normal surface 

tractions which in turn can also simulate any loadings that are applied at dis- 
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tances away from the crack by means of the principle of superposition. 

SHEAR IMPACT 

Suppose that the crack in Figure 1(b) is now sheared suddenly by a pair of 
shear stresses of magnitude such that the upper and lower crack surfaces 
move in the opposite direction. This creates a deformation field that is skew- 
symmetric with respect to the y=0 plane. Following the footstep laid out in the 
previous example on normal impact, the Laplace transform of the transient bound- 
ary conditions on the x-axis inside the crack are 

(t* ) (x,o,p) = - -fp; i<^) (x,o,p) = 0, 0<x<a (34) 

1 H 1 

and the skew- symmetric conditions outside the crack are given by 

(u*) (x,o,p) = 0; (cT*) (x,o,p) = 0, x>a (35) 

X 1 y 

Continuity of the stresses across y=h is expressed by 
(a*) (x,h,p) = (a*) (x,h,p) 

y 1 y 2 

(36) 

(t* ) (x,h,p) = (t* ) (x,h,p) 
xy 1 xy 2 

while the displacements are also required to be continuous, i.e., 

(u*) (x,h.p) = (u*) (x,h,p) 

X 1 X 2 

(37) 

(u*) (x,h,p) = (u*) (x,h,p) 

y 1 y 2 
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Jntzg^aZ A.zpAeJ>e.ntat^yi6. Under the above considerations, the following 
wave potentials <J>^(x,y,p) and i|j^(x,y,p) are selected: 


4>f(x,y,p) = f / [A^^^(s,p)e + A^^^(s,p)e^^^^] sin(sx)ds 

(38) 

'i»f(x,y,p) =|J [B^^^(s,p)e + B^^^(s,p)e^^^^] cos(sx)ds 


for the cracked layer and 


o °° /-I \ "Y-1 oy 

4>f(x,y,p) = f / C^'^(s,p)e sin(sx)ds 
0 

(39) 

4'|(x.y,p) = I / C^^^(s,p)e cos(sx)ds 


for the outside material. 

Equations (38) and (39) may now be substituted into the Laplace transform 
of the stresses and displacements in equations (1) and (2). Making use of the 
conditions in equations (34) to (37), the solution can be expressed in terms of 
the functions A^^^, which are related to a single unknown B(s,p) 

as shown by equations (III.l) in Appendix III. The function B(s,p) is governed 
by the system of dual integral equations 


/ B(s,p) cos(sx)ds = 0, x^a 
0 


( 40 ) 


7TT, 


/ sFjj(s,p) B(s,p) cos(sx)ds = 


x<a 
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The function FjjCs»p) is related to Fj(s,p) in equation (22) as 


( 41 ) 


where 


A =_£ i _. [ 6 ( 2 ) + 6 ^ 3 ) “"^^11 

^II 2^Y21‘-P ® 


-2(Yn+Y2i)h 


. 3O) -’Til"] 


“2Ymh 


The other parameters such as ic^, Aj, 
fined earlier for the case of normal 


g(l)^ g(2)^ etc., are the same as those de- 
impact. 


A solution to equations (40) is again found by application of the Copson's 
method [5]: 






(43) 


provided that 4>|j(C>p) satisfies a Fredholm integral equation of the second kind: 


1 

$fl(?.p) + / ^fi(n,p) Kjj(?,n,p)dn = 


(44) 


whose kernel Kjj(^,n,p) takes the form 

00 

Kii(?,n,p) = P) " (^^) 

Mode II dynamic 4^e44 intmi>ity iacjtofi. As in the case of Mode I, the 
asymptotic expressions of the dynamic stresses in the Laplace transform plane 
are first obtained in terms of r-| and 9^ defined in Figure 2. The results are 
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io*) (r,,e,,p) = - — r sin ^ (2 + cos ^ cos -j-) + 0(rp 

X ^ I I ^ ^ ^ ^ I 

koCP) 

(<^v) i^Tyh»P) = sin j- cos ^ cos -^+ 0(rp 

y 1 I I ^ ^ c ^ 1 

k|(p) 01 , , 

i<y*) (r, ,e, ,p) = - 2v, sin ^ + 0(rp 

z ^ 1 1 ^ 1 ^ I 

kp(p) 01 ®i 

(t* ) (r^,0^,p) = -z— sin cos j- cos -g- + 0(r|) 


( 46 ) 


/2r^ 

with k${p) being the only quantity that depends on time through the parameter p: 


k|(P) = X, 




(47) 


Equation (10) is then applied to invert the Laplace transform of the stress in- 
tensity factor in equation (47). This gives 


k2(t) = 


Tq»^ r 4>fl(l>P) 

HT l^—p 


eP^dp 


(48) 


in which $*j(l»p) is computed numerically from equation (44). 

Figures 10 to 12 display the values of $|j(l>p) as a function of the normal- 
ized quantity C 2 -j/pa for various values of a/h and ^ 2 /^! while “ ^2 ” 

= P 2 are used for all cases. With a knowledge of $^j(l,p)» the integral in 
equation (48) may be evaluated by a procedure outlined in Appendix II. In gen- 
eral, k 2 (t) increases with time reaching a maximum and then decreases to the 
static value for sufficiently large time. The trend is very similar to k-jCt) 
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for the case of normal impact in that a higher value of k 2 Ct) is obtained when 
the modulus of the surrounding material is less than that of the cracked layer, 
i.e., ^ 2/^1 Comparing the results in Figures 6 and 12, it is seen that for 

< 1. normal impact yields a higher crack tip stress intensity factor than 
shear impact, i.e., k^(,t) > k 2 (t). The opposite is observed when ^ 2/^1 ^ i-e., 

k 2 (t) > k-|(t). The curves in Figures 14 and 15 for k 2 (t) also show the absence 
of a small fluctuation for small time which was present in Figures 7 and 8 for 
k^(t). This is because the influence of the reflected incident shear wave from 
the interface is considerably weaker even for the ratio of a/h = 2.0. 

CONCLUSION 

As composite materials are currently being applied to major primary structure 
designs, it is necessary to have an in-depth understanding of the mechanical be- 
havior of these materials, particularly with reference to the allowable applied 
load both statically and dynamically. This investigation is concerned with the 
dynamic stress distribution around a crack . embedded in the matrix of a unidirec- 
tional composite. The time-dependent loading can be of a general nature applied 
in an arbitrarily direction with reference to the crack plane. For those com- 
posites which fail predominantly by matrix cracking under impact, the present re- 
sults can be used effectively for determining the ability of the composite to 
absorb energy and to withstand load prior to total destruction. 

The other modes of failure such as fiber breaking and/or debonding of fibers 
from matrix are not treated but may be significant in other composite systems. 

The redistribution of dynamic stresses in these cases may also be analyzed such 
that their individual contribution can be assessed quantitatively. These cases 
will be left for future investigations. 
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APPENDIX I: EXPRESSIONS FOR AND A^^'^(s,p), — , C^''^(s,p) 

IN NORMAL LOADING 

This section gives the expressions for in equations 

(25) in terms of the variables s, p and the material constants 


= -s[(l - - T# 


'‘2,0^,, ^2rlf22 


’ ''21 y.| '■2 c^2 2 Y 22 


)] 


(2) - cm . % ^ / p^ w ^21~^22 

Ul ^21 y^ 2 c22 


)] 


a 


(3) . i (.2+^2 1 . 1!2 ^ _EL , ^'~^11^22 ,j 

2 ^21' u, t 5^ 's'-y,2V22 


aW - i (s2_2 1 . il2 r,2 ^ (!!!M22)] 

a 2 ^21 ' y^ ^ 2^ ^^"'''l2'*'22 


a 


C5) = . 


1 /.2. 2 X . ^2 r-2 . P^ "^12^2^n 


(I.l) 


a 


(6) = . 1 


Pc 


0 (s''+Y 2 i) y7 


■^12^21 


2c22 s “Y-I2T22 


(7) .r/, ^Z^ ^2 , P^ X/ ^H‘^12 XT 

O s[{l - -)yi, - ^ ( 2 S|j)(s^-V, 2 Y^P 

(8) - _crn . ^ / P^ w ^n'^12 x-1 

“ ■ " " y^^^ll " y^ ^2c|2^ s''-y^2>'22 ^ 


in which y. • is given by equations (16). 

* J 
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The functions are related to the single function A(s,p) 

as follows: 


aO)(,,p) , M - SY„ 


A<^>(s.p) = - M|^ [sv„ .1 (e<D 


-(Yn+Yoi)h 


-2Ynh 


, eC3) 


Bn)(s.p) = bO) a(')(s.p) . b' 2) A<2)(s,p) 

B<2>(s.p) = , g(4) ^<^H-V2l)'> 

n\ 0^12*^ "Ynh /^^ 

^ s^-Yi2Y22 ^^^^‘"^11^22^ ® ^ 

+ (s^+Yi-|Y 22) e A^^^(s,p) + s(y2i-Y 22) e B^^^(s,p) 
- s(Y2i+Y22) 

/o\ a^ 22 -y- - h /-• \ 

‘ s(Tii+Y 22) A^^^(s,p) + (s^-Yt2T2i) e B^^^(s,p) 
+ (s2+Yi2T2i) B^^^(s,p)] 


(1.2) 
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APPENDIX II: METHOD FOR EVALUATING THE DYNAMIC STRESS 

INTENSITY FACTOR EQUATION (31) 

The integral in equation (31) is basically of the form 

The Bromwich path, Br, consists of an infinite line parallel to and to the right 

* 

of the imaginary axis in the complex p-plane. The function f*(l,p) is considered 
to be known for discrete values of p. There are a number of available methods 
for finding g(t) as a process in the Laplace inverse transform. The method 
adopted here can be found in [9,10j. 

The integral f*(l,p)/p in equation (II. 1) is first evaluated at the points 

p = (l+n)6, n = 0,1,2, — (II. 2) 

in which 5 is a real and positive number. According to equations (9) and (10), 
f*(l»p)/P niay be written as 

= / g(t)e’P*dt (II. 3) 

P 0 

The above infinite integral is now transformed to a finite integral on the in- 
terval [-1,1] by making the substitutions 


f*(l,p) stands for $|(l,p) in normal impact and $jj(l,p) in shear impact and 

they are calculated from the Fredholm integral equations of the second kind, 
namely equations (27) and (44). 
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(II.4) 


X 



- 1 


and 

G(x) = g[tCx)] = g[- 1 1og(^)] (II. 5) 

Therefore, equation (1 1. 3) becomes 

■ t*Q^ . lltP . ) i S . ] . = ^ (^ 1 +x)" G(x)dx (II. 6) 

in which G(x) can be expanded in series form consisting of Legendre polynomials 
Pj^(x) which are orthogonal on the interval [-1,1], i.e., 

00 

6(x) = I C. P.(x) (II. 7) 

i=0 ^ ’ 

Similarly, the function (1+x)^ in equation (II. 6) may also be expanded in the 
form 

(1+x)" = I D. P.(x) (II.8) 

i=0 ’ ^ 

such that 


•'i 


2"(2i+l) 


n(n-l)---Cn-(i-l)] 

(n+l)(n+2)---(n+i+l) 


(II. 9) 


Putting equations (1 1. 7) and (1 1. 8) into (1 1. 6) and applying the orthogonality 
conditions for the Legendre polynomials, the following sum is established; 
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(II. 10) 



Thus the coefficients C. may be found with given by 


Cq = f*0.5) 


(II. 11) 


For a finite number of N coefficients, a partial sum for G(x) in (II. 7) is ob- 
tained and an approximate evaluation of gCt) can be made since from equation 
(II. 5) 


g(t) = I C.P. C2e”''^ - 1] 
i=0 ^ ’ 


( 11 . 12 ) 


The parameter 5 is chosen such that g(t) is best described for the range of t 
considered. 
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APPENDIX III: EXPRESSIONS FOR A^^'^(s.p), — , C^^*^(s,p) 

IN SHEAR LOADING 


In the skew-symmetnc problem, the unknown functions in equations (38) and 
(39) can also be expressed in terms of a single unknown B(s,p) in accordance with 
the following relationships: 


.0), 


AV^(s.p) = . 


,( 4 ), 


-2y 


^II 




a(2)(s_p) , 

“ll 


-2y 

[sYgiS 



(sO) . B(3)e'^^21 


h 

) 


+ j Cs^+Y2i ) e 


(Yii+Y2i)h^ 


b('>(s,p) = - a(i)(s.p) - a'^) 


(s,p) 


b(2)(s,p) = - 3(3)e“^^n"^2l)»^ ,(l)(3^pj . ,(4)^^^ir^2l)h ^(2)^^^^^ 

" s"^ 12Y22 '^^"'“^11^22^® A^^^(s,p) 

+ (s^+Y-,-|Y22)e^^^ A^^^(s,p) - s(Y2-|-Y22)e B^^^(s,p) 

s ( y2i"*’Y 22)® ^(s,p)l 
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(III.l) 


a(1)(s,p) 

+ s(Yi2+Y-]-|)e^^^ A^^^(s,p) + (s^-Y2‘|Yi 2)® B^^^s.p) 

+ (s^+Y2iYi 2)® B^^^(s,p)] 

where Ajj is given by equation ( 42 ). 
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cr(t) 



Composite 


Figure 1. Fiber- reinforced 





Model 


unidirectional composite subjected to angle impact 



Figure 2. Stress element near crack in matrix of fiber-reinforced 
composite 
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Figure 5. Variations of with C 2 -|/pa for 
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k|(t)/a-| 


2.4 h- 


/^2/A^I = 0.1 


a/h = 1.0 
v\ - 1^2 “0.29 
P\ =Pz 



Figure 6. Dynamic stress intensity factor k, (t) versus time for 
a/h =1.0 
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Figure 7. Dynamic stress intensity factor k, (t) versus time for 
V‘ 2 / “ 0*1 
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Figure 8. Dynamic stress intensity factor k^(t) versus time for 




Time i 


Figure 9. Applied stress as a general function of time 



^21 ^ P® 

Figure 10. Variations of $jj(l,p) with C2-|/pa 
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a/h = 2.0 


= 0. 1 

V| = 1^2 ® 

P\ =Pz 
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0 2,0 3.0 4.0 5.0 

O21 / PQ 

Figure 11 . Variations of «jj(l»p) with C 2 -|/pa 


a/h = 0.5 


P-z/H ' 10.0 
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Figure 12 . Variations of $jj(l,p) with 



k^(t)/T 



0 2.0 4.0 6.0 8 . 

C 2 1 t/a 


Figure 13. Dynamic stress intensity factor k«(t) versu 
a/h =1.0 
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2.0 


a/h = 2.0 



Figure 14. Dynamic stress intensity factor l< 2 (t) versus time for 
yg/u-j = 0.1 
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k2(t)/r 


1.2 


a/h = 0.5 
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Cg I t/ci 


Figure 15. Dynamic stress intensity factor k^Ct) versus time for 
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COMPUTER PROGRAMS 


Nonmat -impact. 


3 

3 

f 3 

3 
3 
3 

' i 

- 5 

<. : 

♦ 

4 

r ♦ 

20 
22 
23 

< 2*» ^ 

31 

33 

r 34 

35 

» 37 

41 

r 46 

47 

5^ 4 

C. 56 

57 
61 
64 

t 64 

72 
72 

74 

75 , 

104 

106 

110 

111 

114 

123 

' 123 

134 
141 
144 

C. 147 

151 
160 
160 

C. 163 

165 
167 
171 

t. 205 

205 
207 
212 


2 

K1 

K2 

K3 

K4 


COMMON/ AUX/Ht^»Pi<l*P*^2tBMUtX»Y 

LP(1)=0.0 

DTA(1) = 0.0 

READ 2t K1,K2,K3»K4 

data sets^to be evaluated 

POINTS 


MATRIX 


TERMS 


= NO. OF 
= NO. OF 
. = NO. OF 
set up data 

AK=K3 

00 5 N=1.K3 
AN = N 

5 pt(n)=AN/AK 
SET UP INTEGRATION 
M=K3-2 
N=K3-1 

A=1^/(3.*A) 

DO 10 K=2fM.2 
10 0(K)=2.*A 

DO 15 K=1»N*2 
15 D(K)=4.*A 

CALCULATE^NONHOMOGENEOUS 

RHS=1.0 ^ 

DO 22 I = 1»K2 
PRINT 9 
9 FORMAT (IHl) 

READ Sl.EMU 
61 FORMATCFIO. 5) 

DO 999 II*1»'<4 

35 N0NUI)=RH^*SQRT (PT (N) > 
calculate kernel MATRICES 
CALL CONST(I) 

DO 20 N=1»K3 
00 20 M=l»’$3 
IF(M-N) 25.30,30 
F(M.N.I) =F CN.M. I) 

ECmT8,I)=FU(I,PT(M),PT(N)» 
CONTINUE ^ 

PRI NT°6*JPT^L) »MCN(L1 

F0RMAT(5X,F8.4.F15.6) 

CONTINUE 

LP<II + 1 >=N 0 N(K 3 ) 

OTA {II+l)=P 

PUNCH^6l.{DTA(IX) .LPCIX) ,IX=1.19) 

22 CONTINUE 
END 


25 

30 

20 


40 


6 

6 

10 

12 

13 

14 
14 
26 
35 


gggSg 5 ? 2 ul 5 K^ipSifpKE.BHU.X.V 

D‘^L=0. 25^ <3” A) 

IF(OEL)40,45,50 
45 SIMP=0. 0 

return 

r=IPl’,S:ilL?j"z’(I,A*3.-CEL. 
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53 Sl=(0EL/3. )♦ CSA+2.*SB+4.»SC) 

61 IF(Sl.EG.O.O) GO TO 45 

62 K=8 

63 35 Sa=SB+SC 

65 DEL=0.5*0EL 

67 SC=Z(I,A+D£L) 

-75 J=K-1 - 

77 DO 5 N=3,J,2 

100 AN=N 

101 5 SC=SC+Z(I,A+AN*DEL) 

113 S2=<0EL/3,)*(SA>2,*SB+4.*SC) 

122 DIFsABS C(S2-S1) /SI) 

125 ER=0.01 ^ ^ 

127 IF(DIF-ER) 30,25,25 

131 30 SIMP=S2 - 

133 RETURN 

133 25 K=2»K 

134 S1=S2 

136 IF <K-2048) 35,35,40 

140 40 PRINT 42,I,A,E 

152 42 F0RMAT(5X,* INT. DOES NOT CONVERGE »,I3,2F9.4) 

152 PRINT 60,X,Y 

162 60 FORHAT( 2F10. 5) 

162 DO 70 J=l,lfl 

166 DIP=J 

167 DIP=DIP/10. 

171 H=Z(I,OIF) 

175 PRINT 60, W 

202 70 CONTINUE 

206 CALL EXIT 

207 END 




SUBROUTINE CHANGE (F, G,D, I) 

7 


REAL F{4,4,1),G(4,4) ,0(4) 
COMMON K1,K2,K3,K4 

7 


7 


DO 10 N=1,K3 

10 


DO 10 M=1,K3 

11 


G(M,N) sF(M,N,I) ♦D(N) 

24 

10 

CONTINUE 

30 


DO 20 N=1,K3 

31 

20 

G(N,N) = G(N,NH-1.0 

40 


RETURN 

41 


END . 

SUBROUTINE LINED ( A ,B ,T , N) 
REAL A(N,N) ,6(N) ,T(N) 

7 


7 


DO 5 I=2,N 

10 

5 

A(I,1)=A(I,1)/A(1,1) 

17 


DO 10 K=2,N 

20 


M=K-1 

22 


DO 15 1=1, N 

23 

15 

T(I)=A(I,K) 

33 


DO 20 J=1,M 

34 


A( J,K)=T(J) 

41 


J1=J+1 

43 


DO 20 I=J1,N 
T(I)=T(I)-A(I,J)»A(J,K) 

44 


55 

20 

CONTINUE 

61 


A(K,K)=T(K) 

65 


IF(K.EQ.N) GO TO 10 

66 


M=K + 1 

70 


DO 25 I=M,N 
A(I,K)=T(I) /A(K,K) 
CONTINUE 

71 

105 

25 

10 

* BACK SUBSTITUTE 

110 


DO 30 1=1, N 

111 

-- 

T(I)=B(I) 

114 


M=I + 1 

116 


IF(M.GT.N) GO TO 30 

121 


DO 30 J=M,N 

122 


B( J) = B(J)-A(J,I)»T(I) 

132 

30 

CONTINUE 

136 


DO 35 1=1, N 
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137 


K=N+1-I 

J.41 


6(K)=T(K)/A(K,K) 

146 


K1=K-1 

150 


IF(K1.£Q,0) GO TO 35 

151 


00 35 Jl=ltKl 

152 


J=K-J1 

154 


T(J)=T( J)-A(J,K)*B(K) 

162 

35 

CONTINUE 

167 


RETURN 

167 

. _ 

END 



FUNCTION BESJO(A) 

3 


IF(A-3. )5,5,10 

5 

5 

a=A»A/9. 

7 


W=l.-2. 2499997*B 

12 


z=B*e 

13 


H=W+1.265620a*Z 

15 


Z=Z*B 

16 


H=H-.3163866*Z 

20 


z=z*a 

22 


W=W+.0444479*Z 

24 


Z=Z»B 

25 


H=W-.0039444*Z 

27 


z=z*e 

31 


8ESJO=W*.00021»Z 

34 


RETURN 

34 

10 

6=3. /A 

36 


W=. 79788456-. 00000077*8 

41 


V=A-.78539816-. 041 66397*6 

44 


Z=B*B 

45 


H=W-. 0055274*2 

47 


V=\/-.00 0C3954*Z 

52 


Z=Z*B 

53 


H=H-. 00009512*2 

55 


V=V+. 00262573*2 

57 


z=z*e 

61 


W=W+. 00137237*2 

63 


V=V-. 00054125*2 

65 


2=2*6 

67 


W=W-. 00072305*2 

71 


V=V-. 00029333*2 

73 


2=z*a 

75 


H=H+. 00014476*2 

77 


V=V+. 00013558*2 

101 


eE3J0=W/SQRTCAl*C0S( V) 

111 


RETURN 

112 


END 


6 


FUNCTION FU(I.A.c) 

COMMON/ AUX/H, P, PKl »PK2 ♦ 5HU, X,Y 

6 


X=A 

7 


Y=d 

10 


IF(A*B)5,10,5 

11 

10 

FUsO.O 

12 


RETURN 

13 

5 

SUM=SIMF(I, 0.0»5. 0) 

20 


ER=0.01 

21 


DEL =5. 0 

23 

20 

UP=OEL+5.0 

25 


AOOL=SIMF(I»OEL,UP) 

32 


DEL =UP 

33 


T£ST=ABS(ADDL/SUM) 

36 


SUH=SUM+ADOL 

37 


IF(TEST-ER) 15,20.20 

41 

15 

FU=SQRT (X*Y)*SUM 

47 


RETURN 

47 


END 
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SUBROUTINE CONST (I) 

3 COMMON/ />UX/H,P,PK1,PK2,BMU,X, Y 

3 PRlsQ.29 

5 PR2=0.29 

6 PK1=SQRT((1.-2.*PR1) / (2.^(1. -PRl))) 

15 PK2=SQRT((1.-2.*PR2) /( Z.^(1.-PR2))) 

24 READ 1,P 

31 1 FORMAT{F10.5) 

31 HH=0.1 

33 HH=10.0- - 

34 HH=5*0 I 

36 HH=4.0 

37 HH=1.0 

41 HH=0.5 

42 HH=2.0 

44 H=1./HH 

45 PRINT 2 ,BMU,FR1,PR2,HH,P 

62 2 FORMAT(/////5X,» MU2/MU1 =*F6.2,* NUl =+F4.2,* NU2 =*F4. 2///5X ,* 

1/H =+F4.2f* C21/PA =*FA.2) 

62 RETURN 

63 END 


5 

5 

6 

10 

11 

13 

22 

31 

40 

47 

53 

54 
57 
64 
71 
77 

10 5 
1.13 
L21 
126 
133 
136 
1.40 
145 
150 
154 
157 
162 
1.65 
170 
173 
176 
202 
205 
2D7 
211 
212 
214 
216 
22 0 
222 
223 
22 6 
2 32 
2.36 
24 2 
251 
255 
263 
26 6 
275 


FUNCTION Z(I,S) 

COMMON/ AUX/H,P,PK1,PK2,BMU,X,Y 

pp=p*p 

C1=PK1»PK1 

C2=PK2^PK2 

GA=s6rT(S*S+C1/PP) 

G3=SQRT(S»S+1./PP) 

GC=SQRT (S»S+C2/BMU/PP) 
G0=SQRT(S*S+1./BMU/PP) 

AA=S*S+l./PP/2. 

AB^l.-BMU 
AC=S*S-GC*G0 
AD=(GB-G0)/AC/PP/2.*BMU 
AE= (G3+GD) /AC/PP/B.’^BMU 
AF= (S^S-6A*GD) /AC/PP/2. *EMU 
AG= (S»S+GA-*GD)/AC/PP/2.+EMU 
AH= (S*S~Ge*GC)/AC/PF/2,*3HU 
AI= (S»S+G8*GC)/AC/PP/2.*EMU 
AJ= (GA-GC) /AC/PP/2.»BMU 
AK= (GAfGCl /AC/PP/Z.*BMV 
Al=-(AB*GB-AO) 

A 2 =AB*Ge-AE 

A3=AA-BMU*S»S-AF 

A4=AA-BMU*S*S-AG 

A5=-AA+BMU*S»S+AH 

A6=-AA+BMU*S*S+AI 

A7=S»CAE*GA-AJ) 

A8 = -S*{AE»GA-At<) 

BA=A1*A6-A2*A5 
EB=A3»A6-S^A2*A7 
SC=A4»A6-S»A2*A8 , 

BD=S*A1*A7-A3*A5 

B£=S*A1*A8-A4*A5 

B1=88/3A 

B2=BC/BA 

B3=B0/BA 

B4=BE/BA 

EA=2.*GA»H 

EB=2.*GE^H 

£C=(EA+£6l/2. 

£0=2. ♦EC 
E1=EXP(-EA) 

E2=EXP(-EB) 

£3=EXP(-EC) 

£4=EXP(-£D) 

DL=B2+B3*E4+E4»£2+B1^£1 
D1=2.*PP/CC/GA/0L 
D2= AA»AA-S»S»GA^GB 
D3=82-B3*£4 

D4=2.*AA+(GB*(31»B4-B2»E2) -S»S*GA)»E3 
05= (AA*AA*S*S*GA*GE) ♦ (34*E2-31*£1> 

- 46 - 



306 
317 
33l 
33 0 


F=01* CD2*O3%0ii+05) 

Z=(F-S) *eESJO(S*X)*EESJO(S*Y) 

RETURN 

END 


5 

5 

5 

5 


C 

C 

c 


16 

16 

1 

30 

30 

2 

34 

34 

40 

99 

44 

44 

101 

65 

102 

65 

67 

- 

73 

73 

300 

112 

10 3 

112 

116 

121 


130 

130 

3 

140 

140 

143 

H»4 

150 

153 

98 

155 

81 

161 

165 

165 

175 

82 

175 

175 

177 

201 

202 

203 

204 

83 

206 

210 

212 

214 

215 

216 

12 

221 

223 

225 

226 
227 
230 

13 

233 

235 

65 

235 

237 

86 

237 

237 

241 

244 

245 

246 

87 


OF JACOei POLYNOMIALS WHICH REPRESENTS A LAPLACE 
INVERSION-INTEGRAL 

DIMENSION A (50) , GLAM (5 01 .PHI (50) fC(*i,50) 

QIMENSICN 0K (101) ,TT (101) 

C0MM0N/2/TI,TF,DT,MN,EK,TT 
READ l.NN.MN.MM 
FORMAT(3I2) 

READ 2,TI,TF,DT 
FORMAT(3F10,5) 

PRINT 99 
FORMAT(lHl) 

CALL SPLICE (GLAM, PHI, MM, C) 

F0RMAT(/////5X.* GLAM PHI *) 

PRINT 102, (GLAM(I),PHI(I) ,1=1, MM) 
FORMAT(5X,F10.5,5X,F10.5) 

M11=MM-1 - 

F0RMAT(////5X,* C(1,J) C(2,J) C (3 , J) 

^^RINT^OS,! (C(I,J),Isl,V) ,J=1,M11) 
FORMAT(5X,F10.5,5X,F10.5,5X,F10 .5,5X,F10.5) 

PRINT 99 
00 10 I=1,NN 
READ 3, SET, DEL 
FORMAT(2F10.5) 

format! /////^ xf ♦BETA s*F5.3,* DELTA =-*F5.3) 

DO 11 L=1,MN 
AL = L 

S=1./(AL+3£T)/0£L , « ^ 

CALL SPLINE (GLAM, FHI, MM, C,S,G) 

F=G*S 

IF(AL-2.)61,82,83 
A(l)=(l.+B£T)*OEL*F 
GO TO 11 

A(2) = ((2. + BET)*OEL*F-A(l))»(3. + a£T) 

GO TO 11 
CONTINUE 
TOP=l. 

L1=L-1 

AL1=L1 

DO 12 J=1,L1 
AJ=J 

TOP=AJ»TCP 

CONTINUE 

L2=2»L-1 

EOT=l. 

DO 13 J=L,L2 
AJsJ 

aOT = (AJ+BET) *BOT 

CONTINUE 

MUL=80T/T0P 

SUM=0.0 

DO lA N=1,L1 

AN=N 

IF(AN-2.) 85,86,87 
T0D=1. 

GO TO 88 
T0C=AL1 
GO TO 88 
CONTINUE 
TOD=l. 

ICH=L1- (N-2) 

DO 15 J=ICH,L1 
AJ = J 

TOD=AJ»TGO -47- 


C(t 



250 

252 

252 

25 ^» 

256 

260 

261 

264 

266 

270 

273 

275 

301 

304 

306 

307 
313 
315 

320 

321 
325 

325 

326 


15 CONTINUE 
88 CONTINUE 

BOC=l. 

JAsLl+N 

DO 16 J=L,JA 

AJsJ 

- --e00=B00*(AJ+B£T) 

16 CONTINUE 
CO=TOD/eOD 
SUM=SUH+CO»A(N) 

14 CONTINUE 

A(L»=MUL*(OEL»F-SUM) 

11 CONTINUE 

CALL JACSERfOEL,A,BET) 
CALL NAMPLT 

CALL 

CALL QIKPLT(TT»BK*101) 
CALL ENDPLT 
10 CONTINUE 
999 CONTINUE 
RETURN 
END 


6 

6 

6 

6 

7 

10 

11 

12 

14 

24 

26 

32 

33 

35 

36 
43 
45 
55 
55 
61 
61 
65 

105 

105 

111 

113 

115 

117 

121 

122 


12 


10 

97 

96 


95 

11 


SUBROUTINE JACSER (D, C , B) 

DIMENSION C(50) ,SF (50) »P(50) 

DIMENSION BK(lOl) ,TT(101) 

COMMON/ 2/TI,TF,OT,MN,EK,TT 

TT(1)=0.0 

BK(1)=0.0 

LM = 1 

T=TI 

T=T+DT 

X=2.*EXP(-0»T)-1. 

CALL JACOBI (MNtXtBfP) 
SF(1)=C(1)*P(1) 

DO 10 L=2,MN 

L1=L-1 

AL=L 

SF(L)=SF(L1) ♦C{L)*F(L) 
CONTINUE 
PRINT 97,T,X 
F0RMAT(/////5X,^ 

PRINT 96 
FORMAT(///5X,» I 
DO 11 1=1,6 


T =*F6.3,* X =»F10.5) 


C(I) 


,5X, 


N 


F(T) 


PRINT 95,I,C(I) ,I,SF(I) 
FORMAT(5X,I2,F10.2,5 


CONTINUE 
LM=LM+1 
BK(LM)=SF(5) 

TT(LM)=T 

IF(T.LE.TF) 60 TO 12 

RETURN 

END 


X,I2,F10.5) 


SUBROUTINE JACOBI (N, X,B ,F8) 

C THIS PROGRAM CALCULATES JACOBI POLYNOMIALS OF 

C - K-1 WITH ARC X AND PARAMETER B 6T -1 
7 DIMENSION PB(N) 

7 AN=N 

10 IF(AN-2.) 1, 2,3 

12 1 PB(1)=1. 

14 RETURN 

14 2 PB(1)=1. 

16 P3(2)=X-e»(l.-X) /2. 

21 RETURN 

22 3 BSQ=B*B 

23 BONE=B+l. 

25 PB(1)=1. 

26 P6(2)=X-B*(l.-X)/2. 

31 00 4 K=3,N 

33 AK=K 


ORDER 
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34 

36 

40 

42 

43 
46 
51 
56 
64 
71 

102 

103 


AK1=AK-1. 

AK2=AK-2. 

K1=K-1 

K2sK-2 

CO l=((2.*AKi)+B)*X 
C01 = ( C2.*AK2)+B)*C01 
C 01 =<( 2 .»AK 2 )+ 80 N£)» (C01-eSO> 
C02=2.»AK2»{AK2 + B)»( (2.*AK1) +B) 
COs2.*AKl*(AKl+8)*((2.^AK2) +a> 
PB(K)=(C01*P8(K1) -C02*PB (K2) )/C0 
RETURN — 

END 


11 

11 

13 

14 
1-1 

15 
21 
23 
23 
23 

25 

26 
26 
32 

34 

35 
35 
41 
43 
43 
43 
54 
65 
65 
65 
70 
72 
72 

74 

75 
77 

100 

100 

102 

103 

103 

106 

106 

111 

113 

114 

115 
121 
121 
123 


10 

11 

12 

13 


14 


15 


16 

17 


18 

19 


1 

101 


SUBROUTINE SPLINE (X» Y ,M »C » XINT » YINT ) 

DIMENSION X(50) « Y(50) »C(4,50) 

IF(XINT-Xa) >1,10,11 
YINT=Y(1) 

RETURN 

CONTINUE 

IF(X(MI-XINT)1,12,13 

YINT=Y(M) 

RETURN 

CONTINUE 

K=M/2 

N=M 

CONTINUE , 

IF(X(K) -XINT)3,14,5 
YI NT=Y(K) 

RETURN 

CONTINUE , , 

IF(XINT-X(K+1))4,15,7 

YINT=Y(K+1) 

RETURN 

YINtJ (x1k 4-1>-XINT)*(CU,K)* (X(K+1) -XINT)^*2+C(3,t<n 
YINT=YIN7-MXINT-X<K) ) ♦ ( C ( 2, K) XINT-X (K) ) »*2<-C ( 4, K> > 
RETURN 
CONTINUE 

IF(X (K-I)-XINT) 6, 16, 17 

K=K-1 

GO TO 4 

YINT=Y(K-1) 

RETURN 
N=K 
K=K/2 
GO TO 2 
LL = K 

K=(N + K) /2 
CONTINUE 

IF(X(K)-XINT)3,14,18 

CONTINUE 

IF(X(K-1)“XINT)6,16, 19 
N=K 

K=CLL+K)/2 
GO TO 8 
PRINT 101 

FORMAK* OUT OF RANGE FOR INTERPOLATION *) 

STOP 

END 


7 

7 

7 

11 

12 

15 

20 

26 

27 

34 

37 

41 


SUBROUTINE SPLICE (X, Y , N ,CJ 

DIMENSION XC50) ,Y(50) , D (50) , P (50 ) , E (5 0) 

DIMENSION A(50,3> ,6(50) ,Z(50) 

MMsM-1 
DO 2 K=1,HM 
D(K)=X(K4l) -X(K) 

P(0 = D(K) /6. 

E(K) = (YCK+1) -Y(K) )/D(K) 

DO 3 K=2,MM 
e(K)=E(K)-E (K- 1 ) 

A(l,2)=-1.-0(1) /0(2) 

A(1,3) = D (1) /D(2) 

A(2,3) = F(2)-P(1) ’A (1,3) 


>C(4,50) 
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r 




44 

50 

51 

53 

54 
61 

_ 65 
70 
7t* 
76 
101 
105 
112 
114 
116 
117 
120 
127 
133 
135 
140 
143 
146 
154 
165 
165 


A(2,2) = 2.*(PC1)+F(2) )-PCl)*A(l,2) 
A(2»3)=A(2,3i/A(2,2) 

B(2)=BC2)/A (2,21 
DO 4 K=3,MM 

A{K,2)=2.*(P{K-1)^P(K))-F(K-1)*A(K-1,3) 

B(K)=B(K)-P(K-1)*B(K-1) 

A(K,3> = P(K)/A(K,2) 

4 B(K) =B( K) /A (K,2) 

0=0 (M-2)/0 (M-l) 

A(M,l) = l,<'Q + A(M-2,3) 

A(M,2)=-Q-A(M,1)»A(M-1,3) 

8(M)=8CM-2)-A(H,l)*B(M-l) 
Z(M)-B(N)/A(M,2) . 

MN=M-2 
DO 6 1=1, MN 

6 Z(K)=8(K)-A(K,3)*Z(K+1) 

Z( 1)=-A(1,2)»Z(2)-A(1,3)*Z(3) 

DO 7 K=1,HM 
Q=1./(6.*D (K)) 

C(1,K)=Z(K)*Q 
C(2,K)=Z(K+1)*Q 
C(3,K)=Y(K) /0(K)-Z(K)^P(K) 

7 C(4,K>=Y(K+1)/0(K)-ZCK+1)*P(K) 

RETURN 

END 
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SheaA ^pact 



3 

3 

3 

3 

3 

3 

3 

4 

5 

20 


20 

22 

23 

24 

31 

33 

34 

35 
37 
41 

46 

47 
54 

56 

57 
61 
64 
64 
72 
72 
7 *+ 
75 

1C4 

106 

110 

ill 

114 

123 

123 

134 

141 

144 

147 

151 

160 

160 

163 

165 

167 

171 

205 

235 

207 

212 


PROGRAM BETA { INPUT , OUTPUT, PUNCHt PLOT , TAPE 99=PL0T) 
real non (LitFCUtit,!) ,G(4,4),0{‘+) ,PT(4) 

REAL 8(4 ) ,C (4) 

REAL LP(50) ,0TA(5Q ) 

£QUIVA;.ENCE (NON, 3) 

COMMON K1,K2,<3, K4 
COMM ON/ A UX/H,P, PK 1 ,PK2,BMU, X ,Y 
LP(l)=C.O 
DTA(1) = C.0 
READ 2,K1,K2,K3,K4 
■ 2 FORMAT(I2) 

♦ K1 = ORDER OF SYSTEM OF EOUATION3 

♦ K2 = NO. OF DISTINCT KERNELS 

♦ K3 = NO. OF DATA POINTS ^ ^ ^ 

♦ K-+ = NO. OF OATA SETS TO BE EVALUATED 
» SET UP OATA POINTS 

AK=K3 

DO 5 N=1,K3 
AN=N 

5 PT(N)=AN/AK 

♦ SET UP INTEGRATION MATRIX 

M=K3-2 

N=K3-1 

A=K3 

A=l./(3. *A) 

DO 1C K=2,M,2 
10 DIK)=2.*A 

DO 15 K=1,N,2 
15 3(K)=4.*A 

D(K3)=A 

♦ CALCULATE NONHOMOGENEOUS TERMS 

RHS=1,C 
DO 22 1=1, K2 
PRINT 9 
9 FORMAT(lHl) 

READ 61,BMU 
61 FORMAT(FlG.5) 

DO 999 11=1, K4 
DO 35 N=1,K3 

35 NON(N)=RHS'‘SQRT (PT (N) ) 

♦ calculate kernel MATRICES 

CALL CONST(I) 

DO 2 0 N=1,K3 
DO 2C M=l,i<3 
IF(M-N)25,3G , 30 
25 F (M,N,I)=F(N,M, I) 

GO TO 2G 

30 F (M,N,I» =FU (I,PT( M),PT(N) ) 

20 CONTINUE 

CALL CHANGE (F ,G,0, I) 

CAlL LINEQ(Gti,C, K^) 

DO 40 L=1,K3 
PRINT 6, PT(L) ,NON (L) 

6 FORMAT(5 X,F6.4, F15 .6) 

40 CONTINUE 

LP( IH-1) =NON (<3 ) 

OTA(II+l) =P 
9 99 CONTINUE 

PUNCH 66 , (OTA ( I X ) , LP ( IX) , I X= 1, 19 ) 

66 FORHAT(2F10 .5 ) 

call LAPINV(0TA ,LP) 

22 CONTINUE 
END 





53 


Sl=(DEL/3.)*(3a+2.*‘SB+4.»3CJ 

< 

61 


IF(Sl.£Q.0. 0) 50 TO 45 


62 


K =3 


63 

35 

3B=SB+SC 


65 


QEL=0.5»0£L 

< 

67 


SC=Z (I*A«-OEL) 


75 


J = K-1 


“ 77 


00 5 N=3, J»2 


loa 


AN=M 

c 

101 

5 

S C=SC*Z(I,Af AN*OEL) 


113 


S2=(0EL/3.) ♦(SA4-2,*Sa+4.*SC) 


122 


DIF = ABS( (S2-S1) /SI) 


125 


ER=C.01 


127 


IF(DIF-ER)30, 25,25 


131 

30 

SIMP=S2 


133 


RETURN 


133 

25 

K=2*K 


134 


31=S2 


136 


IF (K-2043 ) 35, 35, ^ 0 


140 

40 

PRINT 42, I, A, a 


152 

42 

F0RMAT(5X,» INT. DOES NOT CONVERGE *,I3,2F9.4) 

{' 

152 


PRINT 60 ,X, Y 


162 

60 

FORMAT(2F10 .5) 


162 


00 70 J=l,10 


166 


0IP = J 


167 


OIP=OIF/10. 


171 


W=Z(I,DIP) 


175 


PRINT 60, W 


202 

70 

CONTINUE 

< 

2 06 


call EXIT 


207 


END 


7 


SUBROUTINE CHANGE ( F,G,D, I) 
REAL F(4,4,l) ,G(4,4) ,0(4) 


7 

7 

ID 

11 


COMMON K1,K2,K3,K4 
00 10 N=1,K3 
DO 10 M=i,K3 
G(M,N) =F(M,N, I) »D(N) 

* 

24 

3C 

1C 

CONTINUE 
CO 20 N=l, K3 


31 
4 0 
41 

2G 

G (N,N)=G (N,N) <-l.G 
RETURN 
E i-0 



7 


SUBROUTINE LINEQ( A,3, T,N) 
REAL A(N,N) ,B{N) ,T (N) 

7 


DO 5 1=2, N 

ID 

5 

A(I,1)=A(I,1)/A(1,1) 

17 


DO 10 K=2,N 

20 


M=K-1 

22 


00 15 1=1, N 

23 

15 

T (I )=A(I ,K) 

33 


00 20 J=1,M 

34 


A ( J, K)=T ( J) 

41 


J1=J + 1 

43 


DO 20 I=Ji,N 

44 


T(I)=T(I)-A{I,J)»A(J,K) 

55 ■ 

20 

CONTINUE 

61 


A (K, K)=T (K) 

65 


IF(K.EQ.N) GO TO IG 

66 


H=K*-1 

70 


DO 25 I=M,N 

71 

25 

A (I,K)=T (I) /A (K,K) 

105 

10 

CONTINUE 


* 3A 

CK SUBSTITUTE 

IID 


DO 30 1 = 1, N 

111 


T (I)=B(I ) 

114 


N=I»-1 

116 


IF(M.GT.N) GO TO 3C 

121 


DO 30 j=M,N 

122 


B ( J) =3( J) -A { J,I ) *T (I) 

132 

30 

CONTINUE 

136 


00 35 1=1, N 



-52- 




137 K=N+1-I 

141 3 (lO=T(<)/A (K,<) 

146 K1=K-1 

15Z IF{K1.£Q.4) GO TO 35 

151 DO- 35 J1 = 1,K1 

152 J=K-J1 

154 T(J)=T(J)-A(J»0»B(K) 

162 35 CONTINUE 

167 RETURN 

16 7 E NO 


3 

5 

7 

12 

13 

15 

16 
20 
22 

24 

25 
27 
31 
34 
34 
36 
ai 

44 

45 
47 

52 

53 
55 

5 7 
61 
63 
65 

6 7 
71 
73 
75 
77 

101 

111 

112 


F'UNCTION 3£SJ0(A) 

IF(A-3.) 5,5,10 
5 3 =A*A/9. 

W=l. -2. 2499997^3 
Z = B'^B 

W=W+1.2656203*Z 

Z=Z*B 

W=M-.3163866*Z 
Z = Z*B 

W=W+ ,0444479*Z 
Z=Z*B 

W=W-.0C39444*Z 
Z =Z*B 

3ESJO=W«-. JC021»Z 
RETURN 
10 B=3./A 

W = .7 9786456-. 00 0 0 0C77*B 
V=A-.78539816-. 04166397*3 
Z=3*3 

W=W-.0G55274*Z 
•>/=•«/-. 0CG03954*Z 

z=z*a 

W = W-.00C09512*Z 
y=^+.00262573*Z 
Z =Z*B 

W=VH-.0C137237*Z 
V = '/-.aC0 54l25*Z 
Z = Z*B 

■W=W- .CCC 723C5*Z 
y = ‘./-.0GC 29333*Z 
Z=Z*B 

W=W+ .00C1-4476*Z 
\/=V*.acC 1355fi*Z 
3£SJ0=W/SQRT(A) *COS(\/) 
RETURN 
END 


FUNCTION FU{I,A,3) 

>^nMMnKi/AMv/u o oiy 4 


6 


COflMON/AUX/H,P, PKl , 

6 


X=A 

7 


Y=B 

1C 


IF(A*3)5,1C ,5 

11 

10 

FU=0-.0 

12 


RETURN 

13 

5 

3UM=SIHP(I,0.0,5.0) 

20 


ER=a.01 

21 


0£L =5.C 

23 

23 

U°=OEL+5 . 0 

25 


AOOL=SIMP(I,DEL,UP ) 

32 


DEL =UP 

33 


T£ST=A0S (AOOL/SUM ) 

36 


SUM=SUN4 AOOL 

37 


IF(TEST-ER) 15,20,20 

41 

15 

FU=SQRT(X*Y)*SUM 

47 


RETURN 

47 


E NO 
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3 

3 

5 

6 
15 
24 
31. 
31 

33 

34 

36 

37 

41 

42 

44 

45 
62 

62 

63 


3U3R0UTINC CONST ( I ) 

COMMON/AUX/H,P, PKltPK2tBHU,X,Y 

P Rl=0.29 

?R2=C.29 

PK1=SQRT ((1.-2.»PR1)/(2.*(1.-PR1))) 
PK2=SQ?T ((l.-2.»PR2)/(2.»(l.-PR2>) ) 
READ 1,F 

1 F0RMAT(F1C.5) 

H H= 0 . 1 
HH=lC.o 

H F=5.0 
HH=4.0 
HH=0 .5 
HH=1.0 
HH=2.0 
H=1./HH 

PRINT 2»BHUtPRl»PR2,HH,P 

2 FORMAT(/////5Xt ♦ MU2/MU1 =»F6.2»* NU 
1/H =*F4.2,* C21/PA =»F4.2» 

RETURN 
E NO 


=»F4.2,* NU2 =*F4.2///5X,* A 


o 


la 

11. 

13 

22 

31 

40 

47 

53; 

54 

57 

64 

7r. 

77 


105 
113 
121 
1 26 

1 33 
136 
140 
145 
150 
154 
15 7 
162 
165 
170 
173 
176 

2 0 2 . 
205 
207 
211 . 
212 
214- 
216 
222 . 


2 22 
223 
2 2 6: 
232 


236 
242 
2 51 . 
2 55 
2 6 3 . 
2 66 - 
2 75 ; 


FUNCTION Z(I,S) 

COMM ON/ A UX/H ,Pt PK 1 »PK2, BMU , X » Y 

pp=p*p 

C1=PK1»FK1 
C2=PK2*PK2 
CC=1 .-Cl 

5A=5QRT(S*S+C1/PP) 

GB=SQRT(S»S+1./PP) 

GC=SQRT<S*5<-C2/3MU/PP) 

GO=SQRT(S*S+l ./BMU/PP) 

AA=S»S+1 ./PP/2. 

AB=1,-BNU 
AC=5»S-GC*G0 
AO= (G3-G0)/AC/?P/ 2 .♦BMU 
AE= (G3+GO)/AC/PP/2 .♦BMU 
uF= (3*S-GA»GO)/AC/PP/2.*BMU 
AG= (S*StGA’‘GD)/AC/FF/2.*aMU 
AH= (S*S-G3*GC)/AC/PP/2.*BMU 
AI= (S*5+G3*GC)/AC/PP/2.»BHU 
A j= (GA-GC)/AC/PP/2.»BMU 
AK= CGA<-GC)/AC/PP/2.»3MU 
Al = - (A3*G3-A0) 

A2=A8*GB-AE 

A3=AA-BNU*S»S-AF 

A4=AA-3NU^S*S-AG 

A5 = -AA<-BMU*S»S«-AH 

A6=-AA+BMU»S*5+AI 

A7=S»(AB»GA-AJ) 

A8=-S*(AB*GA-AK) 

3A=A1*A6-A2*A5 

3B=A3*A6-3*A2*A7 

3C=A4*A6-S*A2*A3 

3D=S*A1* A7-A3»A5 

3E=S*A1*A3-A4*A5 

31=3B/EA 

9 2 = 3 C/ BA 

33=30/3A 

34=BE/BA 

£A=2.*GA*H 

EB=2.*Gb*H 

E C=(EA4-EB)/2. 

ED=2 .*EC 
E1=EXP(-EA) 

E2=EXP(-EB) 

E3=EXP(-EC) 

E4=£XP<-E0) 

DL=B2+B3*E^-Bh*E2-Q1*E1 
01=2 .♦PP/CC/G3/DL 
02=AA»AA-S*S*GA»GB 
03=32-33»E4 

D4=2.»AA*(G3*(31*B4-B2»B3) -S^S’^GA) *E3 
D5 = (AA’^AA + S*S''3A*GE)»(B4»E2-Bi»El) 
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306 

313 

33G 

330 


|l?r4??l?l:8(s8r.-8£SgOlS 




RETURN 
END 


) 


■ 5 
5 
5 
5 
5 

16 

16 

ZZ 

30 

3^» 

34 

4u 

44 

44 

65 

65 

67 

73 

73 
112 
112 
116 
1 21 
130 
130 
140 
140 

143 

144 
13 ■: 
153 
15 5 
161 
165 
165 
175 
175 
175 
177 
201 
202 
223 
204 
206 
210 
212 
2 14 
215 
2 16 
221 
223 

225 

226 
227 
2 30 
233 
235 
235 
237 
237 
237 
241 

? iiU 

245 

246 


C 

C 

C 


INVERSION. INTEGRAL 

OINENSICN A (50) .GLAM (50 ), PHI (50 ) .C(4.50) 

k -r r. 4 1 


DIMENSION |K<101) .TT(lCl) 
COMMON/2/TI.Tr ,OT ,MN,3K, 


1 

2 

99 

101 

102 

300 

105 

3 

99 


PHI 


♦) 


=*F5.3t* OELTA =»F5.3) 


3 1 

32 

33 


12 


13 


85 

36 

37 


. TT 

READ l.NNtMN^MM 
F0RMAT(3I2) 

READ Z.TI.TF.OT 
r ORHAKoFlO .5) 

PRINT 99 . 

call'^sfliceiglam.phi.mm.C) 

PRINT IDl 

FORMAT(5X,F10.5.5X.F10.5) 

M11=MM;1 

F0RMAT(7///5X,* C(1*J) C(2.J) 

^PRiNriC3,((C(I.J) .I=1.4).J=i»Mll> . P,- 5, 

PORMAT(5X»F10*3.5XtF10«p.5X*F10.5.5X»Flli«5i 

PRINT 99 
DO 10 1=1.HN 
READ 3,BET»0EL 
FORHAT(2F10.5) 

PRINT 93 .SET.OcL^^.^^ 

FORMAT(/////5Xt ♦8siTA 
DO 11 L=1.MN 
AL=L 

C ALl'^SPLINc (GLAM. PHI t MMt C. S . G) 

P=G*S 

IF(AL-2.)81. 02,83 
A ( 1) = (!• «'3cT ) ♦OtL*F 

A(2l=((2.43ET)»0EL*F-A(l)l»(3.*-BET) 

GO TO 11 
CONTINUE 
TOP=l. 
l 1=L-1 
AL1=L1 
DO 12 Jsl.Ll 

AJ=J 

top=aj»top 

CONTINUE 
L2=2*L-1 
30T=1. 

D C 13 j=L.L2 

30T= (AJ^BET ) ♦BOT 
CONTINUE 
MUL=BOT/TOP 
3UM=C.!: ^ ^ 

DO 14 N=1.L1 
AN=N 

I F( AN-2. ) 85.86.87 
TOO=l. 

goto '88 

T03=AL1 
GO TO 88 
CONTINUE 
TOO=l. . 

ICH=Ll-(N-2) 

DO 15 J=ICH»L1 
AJsJ 

TO 2= AJ*TOO 


C(3, J) 


C(4 
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250 

15 

GONT INUE 

252 

8a 

CONTINUE 

252 


3 0 0= 1 . 

254 


J A=lH-N 

256 


DO 16 J=L,JA 

260 


A J= J 

261 


30D=BOD* (Aj+BET) 

264 

15 

CONTINUE 

266 


CO=TOD/BOO 

273 


SUM=SUM+CO*A IN) 

273 

14 

CONTINUE 

275 


A (L) =MUL* (0£L*F- 

301 

11 

G ON T INUE 

304 


CALL JACSER(DEL', 

306 


CALL NANPlT 

337 


CALL QIKSET (6.C , 

313 


GALL QKSAX(3,3I 

315 


gall QlKPLT(Tr,3 

320 


GALL ENDPLT 

321 

ID 

CONTINUE 

325 

999 

CONTINUE 

325 


RETURN 

3 26 


END 


6 

6 

6 

6 

7 

Ij 

11 

12 
14 
24 
26 

32 

33 
35 
3 o 
45 


45 
5 5 
55 
61 
61 


65 
ia5 
105 
111 
113 
115 
117 
121 
1 22 


SUBROUTINE JAGS ER ( 0, C ,B ) 

DI'IENSION C(50) tSF(5Q)tP(50) 

DI.-IENSICN BK(lGl) ,TT(101> 

G OHM ON/2 /TI ,TF, CT ,KN,3K, TT 

T T{l)=C.O 

BKtl)=C.O 

LM=1 

T=TI 

12 T=T4-DT 

X=2.*£XP (-□♦T)-l. 

GAlL JACOBI (MN, X, 3, P) 

3F(1)=C(1I*P(1) 

DO IG L=2,MN 
Li=L-l 
A L = l. 

SF(l)=5F (l1)<-G(L) ♦?(!> 
lu GOUriNUt 

=RINT 97»T,X 

97 FORMAT(/////5X, * T =*F6.3,* X =*F1C.5) 

=RINT 96 

96 FORMAT(///5X,» I C(I) ♦♦SX,* N F(T) »> 
uO 11 1=1,6 

DRINT 95 ,I,C(I) ,I,SF(I> 

95 FORMAT(5X,I2,F10. 2,5X,I2,F1C .5) 

11 CONTINUE 
LM=Lh+l 
3K{LM)=SF(5) 

TT (lM)=T 

IF(T.Lt.TF) GO TO 12 

RETURN 

END 


7 

7 

1C 

12 

14 

14 

16 

:i 

;2 

23 

25 

2b 

21 

23 


SUBROUTINE JACOBI (N,X,B,P3) 

C THIS PFOGRAM CALCULATES JACOBI POLYNOMIALS OF ORDER 

C <-l WITH ARG X AND PARAMETER B GT -1 

DIMENSION PB(N) 

AM=N 

IF ( AN-2. ) 1,2,3 

1 PB(1)=1. 

RETURN 

2 P3(l)=l. 

PB{2)=X-a*(l.-X)/2. 

RETURN 

3 3S0=B»B 
30NE=B+1 . 

PB{1>=1. 

=3(2)=X-3»(l.-X)/2. 

DO 4 K=3,N 
AK=K 
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3^+ 
36 
40 
■ 42 
43 
46 
51 
56 
64 
71 
102 
1Q3 


aki=ak-i. 

AK2=AK-2 . 

K1=K-1 

K2=K-2 

C01=((2.*AKl)v3)*X 

3 01= { !i* *Al<l!tioNE)^(C01-aSQ) 

^02 = 2»*AK2*^(AK2«'3) ♦( (E^’^AKl) ♦•3) 
3oii:^AKlMAKi;a)M(2.»AK^^ 

P3(K)=(C01*PB(K1) -C02*PB(K2» )/C0 

RETURN 

ENO 


11 

11 

13 

14 

15 
15 
21 
23 
23 
23 

25 

26 
26 
32 

34 

35 
35 
41 
43 
43 
43 

54 
65 

55 
65 
7C 
72 
72 

74 

75 
77 

loa 

IOC 

1G2 

103 

103 

106 

106 

111 

113 

114 

115 
121 
121 
123 


13 

11 

12 

13 


14 


15 


16 

17 


13 

19 


101 


SU3R0UTINE SPt-INE (X,Y »Mt C»XlNT ♦YINT) 
DIMENSION X(53) ,Y(53),C(4,50) 

IF(XINT-X (1) ) It 10 ♦ 11 
YINT=Y(1) 

RETURN 

CONTINUE 

IF(X (M) - XINT) It 12 1 13 

YINT=Y(D 

RETURN 

CONTINUE 

<=M/2 

N=M 

CONTINUE , ^ _ 

IF(X (K)-XINT) 3tl4t5 
YINT=Y(K ) 

RETURN ^ , 

CONTINUE _ _ 

IF(XINT-X(K<-l))4t I5t7 
YINT = Y(K + 1) 

RETURN 

RETURN ^ 

CONTINUE 

I F ( X ( K-1 )-XlNT)6tlctl7 

KsK-1 

GO TO 4 

YINT=Y(K-1> 

RETURN 
N=K 
K=K/2 
GO TO 2 
LL=K 

<s(N+K)/<; 

CONTINUE , . 

IF(X {K)-XINT) 3tl*t,18 

5P'l5^r-i.-XINT.6. 16.19 
N=K 

K=(LL+K)/2 
GO TO 8 

?ORMAT(“oUT OF RANGE FOR INTERPOLATION M 

STOP 

ENO 


7 

7 

7 

11 

12 

15 

20 

26 

27 

34 

37 


DIMENSION A(50t3)tB(50)tZ(50) 
MM=N-1 

a°t<?=xlK+M-x{K) 

^ (k 1 =*(y'(K ^1 » - YCK) ) / 0 (K) 

DO 3 K= 2 tMM 
3 (<)=£(K)-£(K-1) 

A (l»2)=-l.-0(ll /O (2) 
A(lt3)=0(l)/0{2) 

A (?,3)=P{2)-Ptl>*A(lt3) 
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kU 

50 

51 

53 

54 
61 
65 
70 
74 
76 

101 
ID 5' 
112 
114 
116 
117 
1 20 
1 27 
133 
135 
140 
143 
146 
154 
165 
165 


tt(2»2) = 2.*(PllM-P(2))-Pm^A(l,2J 
a (2» 3»=A (2t 3) /A (2 ♦ 2) 

3 (2)=0C2)/A(2,2) 

DO 4 K=5,iH 

A CKt2)=2.*(P(K-l> +P(K))-P(K-1) *A (K-1, 
B {K)=B(K)-P(K-1)*B(K-1) 

A (K, 3)=P(K)/A{<,2) 

3 (K)=3(K)/A{K,2) 

Q=0 (M-2) /Q(M-l) 

A (M, l)=1.4Q+AlM-2,3) 

A (M,2) = -Q-A(M,1)*A(M-1»3) 

3 CM)=3(H-2)-A(M,l )*3{«-l) 

Z (M)=a(M)/A (Mt2) 

MN=M-2 
DO 6 1=1 fMN 
K ** M* I 

zTk) =3(K) -A (K»3) *Z (K>1) 

Z (1) =-A(l,2)*Z(2) -A(lt3)*Z(3) 

DO 7 K=lt^M 
Q=l./(6. ♦D(K)) 

3(1,K)=Z (K)»Q 

C (2tK)=Z (K+1) *0 

3 (3,K) = Y( K) /D(K) - Z (K) *P(K) 

C I4» K)=Y (K*-l» /O (K) -Z (K>1)*P (K) 


RETURN 
E NO 


3 ) 
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